home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power Programmierung
/
Power-Programmierung (Tewi)(1994).iso
/
magazine
/
drdobbs
/
1991
/
05
/
pseudo.asc
< prev
next >
Wrap
Text File
|
1991-03-15
|
5KB
|
157 lines
_A FAST PSEUDO RANDOM NUMBER GENERATOR_
by W.L. Maier
[LISTING ONE]
/******************************************************************************
* Module: r250.cpp Description: implements R250 random number generator,
* from S. Kirkpatrick and E. Stoll, Journal of Computational Physics, 40,
* p. 517 (1981). Written by: W. L. Maier
******************************************************************************/
#include <stdlib.h>
/**** Static variables ****/
static unsigned int r250_buffer[250];
static int r250_index;
/**** Function prototypes ****/
void r250_init(int seed);
unsigned int r250();
double dr250();
/**** Function: r250_init Description: initializes r250 random number
generator. ****/
void r250_init(int seed)
{
/*---------------------------------------------------------------------------*/
int j, k;
unsigned int mask;
unsigned int msb;
/*---------------------------------------------------------------------------*/
srand(seed);
r250_index = 0;
for (j = 0; j < 250; j++) /* Fill the r250 buffer with 15-bit values */
r250_buffer[j] = rand();
for (j = 0; j < 250; j++) /* Set some of the MS bits to 1 */
if (rand() > 16384)
r250_buffer[j] |= 0x8000;
msb = 0x8000; /* To turn on the diagonal bit */
mask = 0xffff; /* To turn off the leftmost bits */
for (j = 0; j < 16; j++)
{
k = 11 * j + 3; /* Select a word to operate on */
r250_buffer[k] &= mask; /* Turn off bits left of the diagonal */
r250_buffer[k] |= msb; /* Turn on the diagonal bit */
mask >>= 1;
msb >>= 1;
}
}
/**** Function: r250 Description: returns a random unsigned integer. ****/
unsigned int r250()
{
/*---------------------------------------------------------------------------*/
register int j;
register unsigned int new_rand;
/*---------------------------------------------------------------------------*/
if (r250_index >= 147)
j = r250_index - 147; /* Wrap pointer around */
else
j = r250_index + 103;
new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
r250_buffer[r250_index] = new_rand;
if (r250_index >= 249) /* Increment pointer for next time */
r250_index = 0;
else
r250_index++;
return new_rand;
}
/**** Function: r250 Description: returns a random double in range 0-1. ****/
double dr250()
{
/*---------------------------------------------------------------------------*/
register int j;
register unsigned int new_rand;
/*---------------------------------------------------------------------------*/
if (r250_index >= 147)
j = r250_index - 147; /* Wrap pointer around */
else
j = r250_index + 103;
new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
r250_buffer[r250_index] = new_rand;
if (r250_index >= 249) /* Increment pointer for next time */
r250_index = 0;
else
r250_index++;
return new_rand / (double)0xffff; /* Return a number in 0.0 to 1.0 */
}
[LISTING TWO]
/******************************************************************************
* Module: rtest.c Description: tests R250 random number generator by
* placing data in a set of bins.
******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
/**** Constants ****/
#define NMR_RAND 5000
#define MAX_BINS 500
/**** Function prototypes *****/
unsigned int r250();
void r250_init(int seed);
/**** Function: main ****/
void main(int argc, char *argv[])
{
/*---------------------------------------------------------------------------*/
int j, k;
int nmr_bins;
int seed;
int bins[MAX_BINS];
double randm;
double bin_limit[MAX_BINS];
double bin_inc;
/*---------------------------------------------------------------------------*/
if (argc != 3)
{
printf("Usage -- rtest [nmr_bins] [seed]\n");
exit(1);
}
nmr_bins = atoi(argv[1]);
if (nmr_bins > MAX_BINS)
{
printf("Error -- maximum number of bins is %d\n", MAX_BINS);
exit(1);
}
seed = atoi(argv[2]);
r250_init(seed);
bin_inc = 1.0 / nmr_bins;
for (j = 0; j < nmr_bins; j++)
{
bins[j] = 0; // Initialize bins to zero
bin_limit[j] = (j + 1) * bin_inc;
}
bin_limit[nmr_bins-1] = 1.0e7; // Make sure all others are in last bin
for (j = 0; j < NMR_RAND; j++)
{
randm = r250() / (double)0xffff;
for (k = 0; k < nmr_bins; k++)
if (randm < bin_limit[k])
{
(bins[k])++;
break;
}
}
for (j = 0; j < nmr_bins; j++)
printf("%d\n", bins[j]);
}